In [1]:
%env PROJ_LIB=/home/pavel/develop/anaconda2/envs/ds/share/proj/
In [2]:
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
Загрузите агрегированные данные о поездках в мае 2016. Просуммируйте общее количество поездок такси из каждой географической зоны и посчитайте количество ячеек, из которых в мае не было совершено ни одной поездки.¶
In [3]:
grouped = pd.read_csv('../grouped_regions_trips_count.csv', parse_dates=['time'])
grouped.head()
Out[3]:
In [4]:
region_trips_count = grouped.groupby(['region'])['trips_count'].sum()
print 'Число ячеек без поездок: %s' % sum(region_trips_count == 0)
Нарисуйте статическую карту Нью-Йорка. Поставьте на карте точку там, где находится Эмпайр-Стейт-Билдинг.¶
In [5]:
lat_min, lat_max = 40.49612, 40.91553
long_min, long_max = -74.25559, -73.70001
empire_lat, empire_long = 40.748817, -73.985428
# EPSG код взят отсюда: http://spatialreference.org/ref/epsg/?search=new+york+&srtext=Search
NYC_EPSG_CODE = 2261
In [6]:
def get_lat_space(num):
return np.linspace(lat_min, lat_max, num + 1)
def get_long_space(num):
return np.linspace(long_min, long_max, num + 1)
In [7]:
%%time
fig = plt.figure(figsize=(15,15))
m = Basemap(llcrnrlat=lat_min, llcrnrlon=long_min,
urcrnrlat=lat_max, urcrnrlon=long_max,
epsg=NYC_EPSG_CODE)
m.arcgisimage(xpixels=2000, service='NatGeo_World_Map')
empire_x, empire_y = m(empire_long, empire_lat)
m.plot(empire_x, empire_y, 'ro', markersize=10)
plt.text(empire_x - 5000, empire_y - 1000, 'Empire State Building', color='black')
plt.show()
Поверх статической карты Нью-Йорка визуализируйте данные о поездках из каждой ячейки так, чтобы цветовая шкала, в которую вы окрашиваете каждую ячейку, показывала суммарное количество поездок такси из неё.¶
In [7]:
lat_space = get_lat_space(50)
long_space = get_long_space(50)
In [8]:
# scaled_data = np.log(region_trips_count + 1)
scaled_data = region_trips_count + 1
data = scaled_data.values.reshape((50, 50)).T
In [9]:
from matplotlib.colors import LogNorm
def plot_static_map_with_colormesh(lat_space, long_space, mesh_data, cmap='hot', xpixels=2000):
fig = plt.figure(figsize=(17,15))
m = Basemap(llcrnrlat=lat_min, llcrnrlon=long_min,
urcrnrlat=lat_max, urcrnrlon=long_max,
epsg=NYC_EPSG_CODE)
m.arcgisimage(xpixels=xpixels, service='World_Topo_Map')
m.pcolormesh(long_space, lat_space, mesh_data,
norm=LogNorm(vmin=1, vmax=np.max(mesh_data)),
latlon=True, alpha=0.5, cmap=cmap)
empire_x, empire_y = m(empire_long, empire_lat)
m.plot(empire_x, empire_y, 'ro', markersize=10)
plt.text(empire_x - 5000, empire_y - 1000, 'Empire State Building', color='black')
plt.colorbar(orientation='vertical', shrink=0.85, label='Trips Count')
plt.show()
In [76]:
%%time
plot_static_map_with_colormesh(lat_space, long_space, data, xpixels=2000, cmap='YlOrRd')
Для большей наглядности выведем ту же карту, но с разбиением ячеек не 50x50, а 200x200. Данные для этого были посчитаны в ноутбуке первой недели и сохранены в файл.
In [86]:
%%time
grouped_hires = pd.read_csv('../grouped_regions_trips_count_hires.csv')
data_hires = grouped_hires.groupby(['region'])['trips_count'].sum()
In [11]:
data_hires.shape
Out[11]:
In [12]:
scaled_data_hires = data_hires + 1
data_hires = scaled_data_hires.values.reshape((200, 200)).T
In [13]:
%%time
plot_static_map_with_colormesh(get_lat_space(200),
get_long_space(200),
data_hires, cmap='hot')
Вставьте интерактивную карту Нью-Йорка — такую, которую можно прокручивать и увеличивать. Поставьте метку там, где находится статуя свободы.¶
In [16]:
import folium
from IPython.core.display import display, HTML
In [17]:
statue_lat, statue_long = 40.689247, -74.044502
map_center = [np.mean([lat_min, lat_max]), np.mean([long_min, long_max])]
In [18]:
def show_map(folium_map):
display(HTML(folium_map._repr_html_().replace('padding-bottom:60%', 'padding-bottom:85%')))
In [21]:
m = folium.Map(location=map_center,
width=800, height=800,
zoom_start=11)
folium.Marker([statue_lat, statue_long], popup='Statue of Liberty').add_to(m)
show_map(m)
Нарисуйте на интерактивной карте Нью-Йорка ячейки так, чтобы их цвет показывал среднее за месяц количество поездок такси в час из этой зоны.¶
In [127]:
mean_hour_trips_count = grouped.groupby(['region'])['trips_count'].mean()
In [149]:
from matplotlib.colors import Normalize, LogNorm, rgb2hex
from matplotlib.cm import ScalarMappable
def plot_folium_map_with_colormesh(lat_space, long_space, density, cmap='hot', fill_opacity=0.5,
cell_filter_func=lambda x: True):
m = folium.Map(location=map_center,
width=800, height=800,
zoom_start=11)
density = density + 1
norm = LogNorm(vmin=np.min(density),
vmax=np.max(density))
mapper = ScalarMappable(norm=norm, cmap=cmap)
for i, lat_low in enumerate(lat_space[:-1]):
for j, lon_low in enumerate(long_space[:-1]):
if not cell_filter_func(density[i][j]):
continue
color = rgb2hex(mapper.to_rgba(density[i][j])[:3])
folium.Rectangle(bounds=((lat_low, lon_low),
(lat_space[i+1], long_space[j+1])),
color=color,
opacity=0.3,
weight=1,
fill_opacity=fill_opacity,
fill_color=color).add_to(m)
show_map(m)
In [150]:
%%time
plot_folium_map_with_colormesh(get_lat_space(50),
get_long_space(50),
mean_hour_trips_count.values.reshape((50, 50)).T,
cmap='YlOrRd')
Чтобы не выбирать из всех 2500 ячеек вручную, отфильтруйте ячейки, из которых в мае совершается в среднем меньше 5 поездок в час. Посчитайте количество оставшихся. Проверьте на карте, что среди этих ячеек нет таких, из которых поездки на самом деле невозможны.¶
In [131]:
len(mean_hour_trips_count[mean_hour_trips_count >= 5])
Out[131]:
In [151]:
%%time
plot_folium_map_with_colormesh(get_lat_space(50),
get_long_space(50),
mean_hour_trips_count.values.reshape((50, 50)).T,
cell_filter_func=lambda x: x >= 5,
fill_opacity=0.8,
cmap='YlOrRd')
На карте мы видим, что все заполненные ячейки (т.е. со средним количеством поездок не менее 5) расположены на суше, т.е. поездки в них не являются невозможными.